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ABSTRACT 

Methods  for  improving  the  calculation  of  energy  values 
for  atoms  by  accounting  for  the  Inter-particle  correlation 
are  presented.   An  approach  using  the  projection  of  an 
intrinsic  state  wave  function  is  first  studied.   A 
formalism  for  the  use  of  intrinsic  state  density  matrices 
is  then  presented  and  applied.   The  special  properties  of 
an  oscillator  based  density  matrix  are  investigated.   As  a 
result  of  these  considerations  a  special  N-body  coordinate 
transformation  is  applied  to  a  Slater  determinant  to  form 
a  correlated  wave  function  and  its  related  one  and  two 
body  density  matrices.   A  calculation  is  carried  out  to 
show  the  improvement  resulting  from  this  approach. 

A  harmonic  atom  approach  with  the  use  of  unitary 
transformations  to  obtain  a  large  class  of  one-body  density 
matrices  and  guaranteed  realizable  two-body  density 
matrices  is  then  formulated.   The  formalism  is  tested  in 
a  calculation  on  Helium. 
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1.   Introduction. 

The  most  prevalent  approach  to  many  body  problems  Is 
the  Hartree-Pock  method  of  approximating  the  solution  to  a 
given  Hamlltonlan  H  by  an  antlsymmetrlzed  product  wave 
function  tj/   composed  of  single  particle  functions  ^.    which 
must  be  determined.   Application  of  the  variational  principle 


(1.1)      5  I   i/*n  T}/  6x  =   0      subject  to  /  ^*^  dx  =  1 


gives  rise  to  a  set  of  non-linear  integro-differentlal 
equations  for  the  (j).  in  which  each  particle  moves  in  the 
average  field  of  the  other  particles.   For  electrons  in 
atoms,  which  we  will  be  dealing  with,  the  average  coulomb 
Interaction  between  one  electron  and  the  others  is  taken 
into  account  neglecting  the  correlation  between  the 
motion  of  the  electrons.   The  best  results  obtained  in 
this  manner  will  be  in  error  by  an  amount  designated  as 
the  correlation  energy  E  ,   The  solution  is  usually  obtained 
by  an  iteration  procedure  starting  from  a  trial  set  of 
functions  (j).  .   Any  given  trial  solution  of  Hartree  Fock 
form  can  be  obtained  from  any  other  similar  function  by 
a  product  of  one  body  unitary  transformations  U(i)  where 
1  represents  the  1-th  coordinate.   Symbolically, 


(1.2)        i/^  =  \Jf     with   U  =  TT  U(i) 

1=1 


The  trial  set  determines  the  potential  which  in  turn 
determines  the  next  set  of  wave  functions.   If  the  process 
converges  it  is  continued  until  a  self-consistent  field 
is  obtained. 

The  wave  function  obtained  in  this  manner  may  give 
fairly  good  values  for  the  energy  but  poor  approximations 
to  the  true  ground  state  wave  function.   The  reason  for 
the  aforementioned  is  that  in  using  the  variational 
principle  first  order  errors  in  the  trial  function  produce 
only  second  order  errors  in  the  energy.   To  obtain  a 
realistic  wave  function  spatial  correlation  effects  must 
be  considered. 

A  number  of  methods  have  been  employed  to  account 
for  the  correlation  energy  [1],  [3]-   One  major  approach 
is  that  of  configuration  interaction  (CI)  [2].   The  true 
solution  must  satisfy  the  correlation  cusp  condition 


(1.5a)  (IJt)   ^^^  =  1 

where  f    is  averaged  over  the  angle  of  r.  ..   For  any 
truncated  configuration  interaction  function  J^^ 


The  numbers  in  [  ]  are  references 


so  that  no  truncated  '^r,j   can  satisfy  the  correlation 
cusp  condition.   The  present  formulations  will  avoid  this 
difficulty  by  involving  a  non-explicit  and  thus 
non-truncated  configuration  interaction.   Completely 
arbitrary  correlated  trial  wave  functions  are  not  in 
general  use  as  in  many  particle  systems  they  would  be 
too  complicated  to  use  and  too  difficulty  to  integrate. 

A  class  of  approximate  wave  functions  which  includes 
correlation  is  given  by  Eq.  (1.2)  when  f^   is  an  interacting 
solution  of  some  related  problem  of  interacting  particles. 
A  still  larger  class  is 

(1.^)        ^-yF(7)U^(l)  ...  U^(N)^Q  da)(7) 

a  generalized  intrinsic  state  model  of  the  sort  used  in 
nuclear  calculations.   In  fact  this  is  the  most  general 
configuration  interaction  expansion  and  thus,  in  a  sense, 
will  include  our  approaches.   To  avoid  the  problem  of 
multiple  integrations  the  formalisms  will  be  developed 
in  terms  of  one  and  two  body  density  matrices  P,  and  Pp 
which,  however,  present  some  difficulties  in  evaluation 
and  the  assuring  of  their  legitimacy  [4].   These  questions 
will  be  examined. 


2.   Symmetry-Breaking  Intrinsic  States. 

Stationary  states  have  symmetry  properties  imposed 
by  the  nature  of  physical  space  which  may  conceal  the  nature 
of  the  specific  system  in  question.   A  prototype  of  this 
Is  the  spherically  symmetric  ground  state  of  a  rigid 
rotator,  the  rotator  Itself  being  quite  asymmetric. 
Our  point  of  view  regarding  this  can  be  that  all  equally 
likely  directions  appear  superimposed  in  the  ground  state, 
giving  a  projection  onto  the  s-state.   If  the  intrinsic 
anisotropic  object  is  a  particle  system,  even  independent 
particles  in  an  anisotropic  potential,  rotational 
averaging  will  produce  new  two  particle  correlations . 
As  an  example  two  electrons  far  apart  have  a  tendency  to 
be  diametrically  opposite  each  other  with  respect  to  the 
nucleus.   The  intrinsic  state  approach  has  been  a  valuable 
model  technique  in  nuclear  physics,  where  short  range 
forces  and  close  packing  produce  a  somewhat  rigid  structure. 
More  surprisingly,  Auffray  and  Percus  [5]  have  shown  a 
rotational-intrinsic  state  to  be  very  good  for  helium,  by 
making  a  configuration  expansion  of  their  result. 

The  projective  technique  that  Is  used  with  the 
intrinsic  state  to  obtain  spin  and  orbital  angular 
momentum  elgenfunctlons  for  a  rotationally  invariant 
Hamlltonlan  is  as  follows.   Let  x-  denote  space 


coordinates  and  Riq-v  a.  space  rotation  with  Euler  angles 

^,   Q,yt     while  u.  denotes  spin  coordinates  and  Si^y 

the  corresponding  spin  rotation  [6].   Define  x.  =  ^^Qy^±> 

(2.1)     u^  =   ^<t)0X^i  '        ^'  =  ^^QX   ^(l)9X^^---^i'^i---^ 


_  1  -  1 


We  now  require 


(2.2)  ^j^(...x^,u^...)  =^^F((t),Q,;^)^'(...x^,u^...)- 

•  sin  Q   6Q   6^   dX    • 
That  is,  we  must  determine  a  weight  function  F((i),9;^) 
which  upon  averaging  over  space  and  spin  rotations  will 
project  out  a  wave  function  ^j^y.  of  definite  angular 
momentum  J  and  its  z  component  M.   It  is  easily  shown 
that  in  order  to  have 

(2.3)  J^^jT^  =  M^j^  and  J^^jj^  =  J(J+l)^jM 

P(<t'^9j/^)  must  satisfy  precisely  the  equations  for  the 
rigid  rotator  whose  solutions  are  well  known  [6] .   In  the 
simple  two  dimensional  case  F{i;),Q,X)   reduces  to 

,   iM0 

(2.4)  f(9)  =  ^  e    . 

It  is  of  interest  to  determine  for  what  potential 
the  projected  wave  function  serves  as  a  solution.   Of 


course  this  will  depend  upon  the  nature  of  the  Initial 
many-body  potential  and  upon  the  allowable  state 
dependence.   Sample  calculations  will  therefore  now  be 
carried  out  using  two  dimensional  cases. 

Consider  a  two  dimensional  product  function 

u{x^,y-^)vix^,y^)    which  we  distort  to  n{ajyi-^,^jy-^)via2^2'^2^2^ 
which  is  not  an  angular  momentum  eigenfunction .   Now  project 
out  an  angular  momentum  eigenfunction 

(2.5)  ^.   =  ?.f  ^  ~J   d^Q   e  "  R((1)q)^(...x^,...) 

0 

where  R((t)Q)^( .  .  .r .  cos  ({)....)  =  f{...r.    cos  ( ({) . +(t)Q )  .  .  . )  =  ^. 
Taking  a  product  wave  function  implies  non-interaction 
between  the  particles.   Distortion  and  projection  gives 
rise  to  a  non-product  wave  function  with  interaction 
between  the  particles. 

2 
V  is  invariant  under  rotation  in  the  sense  that  if 

(2.6)  v2f(x,y)  =   g(x,y),   V^RCc)))  f  (x,y )  =  R{^Q)g{x,y) 
=   V'^  f(x',y')  =  g(x',y')  . 

So  that  if  the  differential  equation  u  and  v  satisfy 
is  known, 

2 

(2.7)  |j^  V^  R((1)q)  uv  =  R((1)q)(v^-  E^)  uv  . 


The  potential  satisfied  by  9.  now  is 

J 


i  i  A*o  -"'*°  t'l^  '!>  «<*o'- 


=  -  E   -  E   + 


1  r""    -iJ^o  , 


The  last  numerator,  with  its  arguments  written  out,  is 
2"     -iJi 


(2.9)  -k   h^o   ^  ^^^0^    (^^^^l^l'Pl^l^  + 


+  V^(a2X2,P2y2^^  u(a-LX-^,p3_y^)  yicL^x^,^^^^'^     ' 
Equation  (2.8)  says 

(2.10)  Tcf)  +  PVuv  =  (E^  +  E^)<1'   • 

Therefore  to  satisfy  a  new  Hamiltonlan,  H',  we  must  have 

(2.11)  PVuv  =  (V'+C)Puv  . 

For  a  rotationally  invariant  potential  V  =  V  we  obtain 
the  identity.   Setting  C  =  0  would  imply  the  unusual 
requirement  that  this  energy  level  be  the  same  for  H  and  H' 
Eq.  (2.11)  is  generally  not  satisfied  except  in  the  highly 
artificial  sense  that 

PVuv   ^ 

(2.12)  V'  =  p^j^  -  C  . 
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Possibly  V  can  be  considered  a  kernel  for  a  sort  of  integral 
eigenvalue  type  equation. 

We  now  apply  this  technique  to  a  two  dimensional  two 
particle  wave  function  solution  of  a  nonlnteracting  harmonic 

oscillator  central  field  potential.   The  ground  state  for 

2     2 
an  unnormalized  single  particle,  u  =  exp(-a  /2  r.,  ) ,  is 

distorted  to 

2  2  2  2 

CL        2         f ^    -,  \CL        2  a         2    /^    T  \a         2    .    2i 

-  —  r        -{q-D—  y  r   _(q-l)—  r   sin   (t) 

(2.13)    u  =   e  e  =   e 

-(|-  +  |(q-l)|-)r2      -  |(q-l)|-  r^cos   2^ 
=   e  e 

-B,r,      -A,r-,cos   2(j)  _ 

=   e  e 


Similarly 

-B„r„      -A  r   cos   2(t)„ 
(2.14)  v^e^^e^^  ^ 


The  resulting  projected  wave  function  which  is  obtained  is 
given  in  equation  (2.l8).   To  see  the  effect  of  this 
distortion  and  projection  in  modifying  the  Hamiltonian  to 
the  form  which  our  new  wave  function  will  satisfy,  rewrite 
the  Schrodlnger  equation  as 

(2.15)    Vj^  =  £  +  ^  k^p^  +  ^  k^p^  +  ^"^   puv      • 
The  result  obtained  is  given  in  equation  (2.19). 
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The  same  procedure  is  also  carried  out  for  a  coulomb 
potential  and  the  results  are  given  In  equations  (2.20) 
and  (2.21), 

The  following  Is  a  summary  of  results  for  the  various 
cases  based  on  the  ground  state  wave  functions: 

1]    Deformed  two  dimensional  harmonic  oscillator 


2   ^         2  2 
_  a 

2 
e 

[I  [I 


(2.16)   (t),  =  P   e   ^     e        '^ 


_   q+1   2  2 

-IT   '^  P  "■---,  ^«n  ^  ,q-l   2  2, 


where  M-  =  2n.   I^  =  modified  Bessel  function  of  first  kind 
Let  B  =  -(q+i)/4  a^,  A  =  q-l/J|  a^ . 

2  P 

(2.17)      V^-   e"   =  1^   [    4B(l+2n)+2(A2+2B^)p2+^^   + 

P 
I 


+  -J±i   {2A+n+8ABA^}] 


n 

2]    Two  non-interacting  particles  in  deformed  two 
dimensional  harmonic  oscillator  potentials 

(2.18)      i      =   P      e"""!^!   ^-V^°^    ^h   ^-V2   ^-A2P2<^os    2(1>2 


-B   p    -B   p         -111®  /-p ^ ^ 

e  I^(  /Ajp^+A|+2A^pjA2P2COs    2(4.2-(t)^)     ), 


with  2 

1  -1      ^iPl    ^^^  ^^1   ^  •'^2P2   ^-^^   ^^2 

(if)  ~  "p   tan 


A 


^p^    cos    (j)^   +  A^pg    cos    2(j)2 


for  a  =  0. 

2 
(2.19)  Vq-E^  =  ^^+  «i32  -  I  k^Pi  -  I  k2P2  +  |fj;  [^A^p2+  ^A^p^ 

4Iq  fB^CA^pJ+A^p^A^PgCOs  2((t)2-<t)i)+B2(A2P2+A;LP^A2P2COS  2((l)2-(t);j_  )l 


^0 


>4^pJ+A2P2+2A^p^A2P2COS  2{^^-^-^) 


3]   Two  non- interacting  particles  in  a  deformed 
hydrogenic  potential. 

-^iPl+^l^l   -a2P2-^2y2 
(2.20   ^     =P   e   ^^   ^e  ^   "^      ^   ^ 


V^  '^^iPl  "^  ^2P2  "  2b-^p^b2p2  cos  ((i)2-<t)]_)  )  i 
-1      ^1^1  ^-^"^  *^1  ~  ^2P2  ^^^   ^2 


(5)  =  cotan 

b-j^p^  cos  (j)^  -  b2P2  cos  ^^^ 

for  LL  =  0. 

"     "       ^^   r   ?       ^1         ?         ^?         2         2 

(2-21)  Vq-Eq  .  I-  [a^  _  _  ,-  a^  -  _^  +  ^^  +  b^ 


Iq  Ja^(b-^p^-b^b2P2COS  ^^-^^)+   a2(b2p2-'bj_t)2p]_cos  ^^'^i^ 


^lPl"^^2p2"^^lPl^2p2  ^°^  ('t>2-^i) 

We  note  that  this  procedure  gives  rise  to  a  rather  odd 
interaction  term  in  the  form  of  a  quotient  of  bessel 
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functions.   It  appears  that  by  Insisting  on  a  pure 
coordinate  potential  we  have  an  Intense  state  dependence 
and  very  complicated  potentials. 

If  we  allow  momentum  dependence  and  dependence  on 
energy  of  state,  a  systematic  development  of  an  effective 
interaction  is  possible.   We  can  write  the  projection  of 
highest  symmetry  generally  as  an  average  over  a  group  of 
unitary  transformations 

(2.22)  P  =  Tu  dmir) 

where  U^  =  U^(l)  .  .  .U^(N) ,  U^U^ ,  =  u^.^,  ,  J  aa^iy)   =  1,  and 
da3(7-7')  =  dooCr'-r)  =  dcuCr)  for  fixed  7'. 
For  example  for  two  dimensional  rotations 

i  0  L 

(2.23)  Uq  =  e^    ^    ,      dcD(0)  =11  . 

Suppose  we  have  an  anisotropic  Hamiltonian  H  =  T  +  A 
such  that  A  is  the  pure  "anisotropic  part"  of  H.  In  this 
case 

(2.24)  jTu^A  U~^  dcD(7)  =  0 

as  well  as  [T,P]  =  0  ,   which  may  be  realized  by  choosing 

Suppose  further  that  Hcf)  =  E(j),  and  define 

(2.25)  f  =  ?^,      ^'  =?'([)  ,     where  P'  =  1-P . 
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Thus 

(2.26) 

so  that 

(2.27) 

and 

Hence 

(2.28) 

where 


A^+T 


'  ^ 

=    E 

* 

.  ^'  - 

,  *'. 

A^+T 


^'  =  (E-T-A^)  ^  A^^  , 


(A^+A^lE-T-A^l'^A^)^  -  (E-T)^ 


H^  =  E^ 


H  =  T  +  A^  +  AgCE-T-A^)  ■'"A2  , 


and  we  can  write 
(2.29)  H^  = 

Now 


H     0 


0   A-,+Tj 


^ 


0 


■=  E 


0 


(2.30) 


so  that 


P  = 


1   0 
0    0 


(2.31)   A^  =  PAP,   A2  =  PAP',   aJ^P'AP,   A^  =  P'AP', 


and 


(2.52)   H  =  P(T  +  A  +  AP'[E-T-P'AP')  ^P 'A)P  +  P'(T+A)P' 


or  finally 
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(2.35)   H  =  T  +  (PAP+P'  AP  '  )  +  PAP  '  (E-T-P 'AP  '  )  "^P  'AP  . 

We  can  retain  at  most  two  body  terms  by  replacing  the 
reciprocal  operator  by  its  expectation 

(2.34)      -^  =  <([)    I    P  '  (E-T-P 'AP')"^P'  I    ^>   /   ^HF'  ;f^    . 


Then 


1 


(2.35)  H  ~    T   +    (PAP+P 'AP ' )    +  -   PAP 'AP     . 
As   a   typical    evaluation,    note   that 

(2.36)  A  =  XZ  A(l)    ,  U^^  -   ^    , 
so  PAPAP    can  be   replaced  by 

(2.37)  PAPA   -    r^Z  U^A(i)U^,A(  j)    dco(7)    dcD(7')     , 

which  gives 

(2.38)  J^Z    (U^A(i)U+)(U^,A(j)U^,)    daAj)    da)(7') 

=  /ZZ    (U^(i)A(i)U  +  (i))(U^,(j)A(j)U^,(j)    da.(7)    d(D(7') 
Hence   defining 

(2.39)  A^(i)    =   U^(l)A(i)U+(l)    ,  [A^]    =/a^    dcD(7)    , 
we   have 

(2.40)  H   ^    T   +^    [A    (1)]     +)  XI    ([A^(i)A^(j)]-[A^(l)][A^(j)]. 
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In  particular  since  we  have  defined  A  such  that 
[A^(l)]  =  0,  then 

(2.41)  H  -^  T  +^XZ  [A  (1)A  (j)]  . 

As  an  example  In  the  case  of  two  dimensional  space 

rotations  with  A  only  coordinate-dependent, 

^  9L  -  :!  ©L 

(2.112)   Aq(1)  =  e^   ^  A(p^,e^)  e  ^   ^  =  A(p^,e^4^) 


so  that 

(2.45)   [A^(i)A^(j)]  =  ^/a(p^,0^+  0)  A(Pj,O^.+  0)  dO 

In  (0.-0.) 

where  A  Is  the  Fourier  component  of  A,  and  this  is  the 
two  body  interaction  which  is  effective.   For  three 
dimensional  space  rotations  very  similar  considerations 
hold,  while  the  inclusion  of  spin  rotations  raises  to 
two  the  rank  of  the  coefficient  of  the  Fourier  terms, 
which  Is  a  very  helpful  generalization. 
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3-   Intrinsic  State  Density  Matrices. 

When  the  intrinsic  state  is  used  as  a  variational 
ansatz,  it  introduces  extra  integrations  to  those  already- 
needed  to  compute  the  many-body  averages  used  to  determine 
expectation  values.   Griffin  and  Wheeler  [9]  use  the 
intrinsic  state  trial  wave  function 

(3.1)  ^(x)  =  J  ^{x;^)    f(p)  dp  ,  i/*{x)    =J^*{x;a)    f(a)  da 

with  f(a)  to  be  determined  in  the  variational  principle 
for  energy  (Eq.  (1.1)),  where  the  Hamiltonian  contains 
particle  coordinates  only.   In  this  case  all  quantities 
dependent  on  coordinates  are  known  and  the  integrations 
over  these  variables  can,  in  theory,  be  performed.   Only 
integrations  over  a  and  p  will  remain  i.e. 

I    f* (x)    H  ^(x)  (dx) 

(3.2)  E  =  - 


TtI/*{x)    ^(x)  (dx) 

=  J  f*{a)   K(a,p)  f(p)  da  dp  /  f  r*  (a)    I(a,p)  f(p)da  dp 

Here  the  "overlap  integral"  I  and  the  "energy  kernel"  K 
are  abbreviations  for  the  expressions 


(5.3) 


I(a,p)|    p  „  fl  1 

>  =  /  (j)  (x.  ...x..;a)  <   >  (t)(x,  .  .  .x^;p)(dx)  . 
K(a,p)J   J  1    ^    [hJ     1    N 


Now  f  is  varied  to  extremize  (3.2).   In  this  manner,  an 
Integral  equation  is  obtained 
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(3.4)         J[K{a,^)    -    I(a,p)E]  f(p)  dp  =  0  . 

This  equation  as  It  stands  Is  very  difficult  to 
deal  with.   Griffin  and  Wheeler  [9]  present  arguments 
that  for  N  large  and  i/   well  represented  by  a  Slater 
determinant  the  overlap  Integral  I(a,p)  can  be  considered 
sharply  peaked .   They  then  approximate  I  by 

-[s(p-a)^] 
(3-5)         I(a,p)  =  e 

and  represent  the  ratio  of  the  kernels  K/l  by  a  simple 
quadratic  function  of  p  and  a 

(3.6)  K(a,p)  =  |Eg+^[2s-4s2(p-a)2]+-|;^o2(^)^>  -lia,^) 

Their  approach  thus  makes  only  very  limited  use  of  the 
generality  of  Eq.  (3-1) • 

The  use  of  Intrinsic  state  wave  functions  will  here 
be  approached  from  the  point  of  view  of  intrinsic  state 
density  matrices.   The  required  one  and  two  body  density 
matrices  for  the  intrinsic  state  are  to  be  found  in 
terms  of  a  basic  independent  particle  one  body  density 
matrix  In  a  fashion  appropriate  to  large  N  and  with  a 
readily  obtainable  correction  series.   Generalizing  to  a 
projection  onto  other  than  the  identity  representation  of 
a  symmetry  group,  we  want  to  determine  the  expectations 
associated  with  the  N  particle  state 
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n  N  -1/2 

(3.7)  ^(l...N)  =  /H  /  dco(7)  F(7)  Tl   U  (i)(N.')     Det  u,  (i), 

J  £=1      ^  ^ 

where  u,  (^)  are  single  particle  orbltals.   The  N  body 
density  matrix  Tptp     Is  thus 

(5.8)  r(l-..N|l'...N')  =^JjF{y)    F*(7')Det  U^(l)uj^(l)- 

•(Det  \J^{£)vi^{£))*   da3(7)  doD(7')  • 

For  a  general  N  body  wave  function  (f)(l...N)  the  associated 
first  order  density  matrix  Is  defined  as 


(5.9)  r  (l|l')  =  N  J  (t)(l,2,  ...,N)  r(l',2,  ...,N)  (dX^) 

where  (dX, )  means  the  Integration  Is  taken  over  all  variables 
except  1  and  1',  (dX^^^  would  omit  1,1', 2, 2',  etc.   The 
second  order  density  matrix  Is  defined  as 

(3.10)  r(l,2ll'2')  = 


N(N-l) 
2 


J    ^{1,2,3,.  ..,'^)    (i)''(l',2',3,  .  ...N)  (dX^2^ 


In  our  case  let  us  consider  our  first  order  density  matrix 
and  expand  the  determinant  by  the  1,1'  row  and  column 
obtaining 
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(3.11)  r(i]i')  =  TTrrrr/-'-/^^^^  ^'^^'^  ZIu/Dupd) 


p,q 


* 


(U  ,(l')u  (1'))   Det(  >   U  (i)u  (1)(U   (l)u  (D)  ) 
^      ^       a?^p   i?l  ^    ^     ^     ^ 

(dX^)  dcu(7)  da)(7') 

n  j[J^F(7)f''(7')U^(1)U*,(1')  ^  ^p^  l^^J^^  '  ^^pq^^" '^  ^ 

•dcD(7)  da)(7')  , 


where 

(3.12)  Q^q(7',7)  =   Det  (  y' U^(2)u^  (2)  (U^  ,  (2)Up  (2)  )""  d2 

Let 

(3.13)  A^^  =  /^s^l)  U*,(l)  U^(l)  Ugd)  dl  =<Ug|u*,U^|u^> 

and  note  that 

.  Tr  In  A 
(3-14)  Qpq(7'7)  =  cofactor  (Ap^)  =  A'^  det  A  =  A'^e 

To  evaluate  the  density  matrix  we  now  make  the 
approximation  valid  for  large  N,  that  the  N  x  N  submatrix 


of  <(u  lu  ,U  lu  ^  is  unitary  as  well.   Hence 

^   ,        rr  *  Tr,,ln  u*  ,u 

(3.15)  Pdll')  =  njj  F(7)F  (7')  e   N    7'  7  y  (i)u*,(l.) 

•  /ZZ  Up(l)u*(2)U*(2)  YZ   U^,(2)Uq(2)u*(l')  d2  dcD(7)da)(7  '  )  = 
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Tr.an  U*  .U 


=   nJJ   F(7)F    (7')    e      ^  ^      ^   U^(l)    U^  ,  ( 1 ' ) 

•  J  ro(l|2)    U^(2)    U^,(2)     ro(2|l')d2   dao(7)    dcD(7'), 

where 

N  ^ 

(3.16)  Todll')    =  ZZ  u^(l)u^(l') 

r=l     ^  ^ 

is  the  one  body  density  matrix  of  the  original  wave  function. 
In  other  words,  we  conclude  that.  In  operator  form 

(3.17)  ^1  =  ^  T-L  f  1 

where 

n  Tr  (  Pr^n  U^ )  . 

(5.18a)  ri(l|2)  =  /  F(7)  e     ^   ^   U^(l)ro(l|2)U^(2)  60,(7) 

or   In  a  slightly  different  notation 

(3.18b)       <l|ril2>    =  j   P(7)    e  °  >'  <l|U^roU^|2>da3(7). 

For  the  two  body  density  matrix,  ["„,  proceeding  In 
the  same  way  we  have  (In  our  second  notation  to  save 
writing  labels) 


(5.19)<1   2|  Tgl    1'2'>  =  |       YZ  /7f(7)    F*(7')e 


Trj^m  U^ 


1^1, V>2>^l''i2 


Tr  In  U 
e     N  7'    u    (l)u      (1)    U    (2)u^    (2)(U    (l>)u^    ( 1 '  )  )'' (U^(2  '  )u      (2'))* 

f  <X'V7' !%><%•  V7' I  V" 
><u      Iu;u^,|u><u       |u;u^,|u     >]    da3(7)    dco(7  '  )  • 
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0 


Again  identifying    \~q  we   obtain 

[J  {F{7)e     N    ^)(F(7')e  ^    ^  ) 

[<iiu^roujb><3iu^,roU^.ii'><2iu^roU^b'><3'iu^,roU^,i2'> 


* 


<3'|U^,roUy,|l'>]  da^(7)  c3cd(7')  , 
or  in  operator  form 

(3.21)  Pg  =  ^  Tg  ^2 

where 

T   o       Tr(  r^ln  U^) 

(5.22)  <1  2|  rp|l'2'>  =  -^^  /  P(7)  e     °    ^ 

/2     J 

•[<i|Uyrou5;ii'X2iu^roU*i2'>  - 


<l|U^ro  U^|2'><2|U^roU*|l'>]  da)(7) 


We  thus  have  a  formalism  for  dealing  with  any  case 
where  we  have  reason  to  believe  that  the  states  of  a 
system  may  be  given  by  projection  of  an  intrinsic  state. 
Given  the  one  body  density  matrix  for  an  intrinsic  state 
we  can  determine  the  correlated  one  and  two  body  density 
matrices  for  the  projected  state.   We  shall  subsequently 
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generalize  this  formalism  to  Include  an  additional 
transformation  and  apply  the  results  to  atoms.   Now  we 
will  consider  the  case  of  a  periodically  bounded  particle 
along  the  x  axis  and  perform  a  two  dimensional  rotational 
average  to  illustrate  the  effects  obtained.   Our  single 
particle  wave  functions  are 

ikx 

(3.23)  (t>i^(x,y)  =  e    5(y)  , 

so  that  our  original  one  body  density  matrix  is 

(3.24)  ro(l|2)  =  ZZ  ^k^^^l^l^  '^l^''2^2^ 


sin  ^  (x^-x^) 
sin  2  (x^-x^) 


for  N  odd  . 

Under  rotation  this  becomes 
N 


(3.25) 


sin  I  (x^cos  0+y^sin  Q-Xgcos  O-ygiin  0) 
sin  "I  (x^cos  0+y^sin  0^x2cos  O-y^sin  0) 


5(y  cos  0-  x^sln  0)  bCygCOs  0-X2Sin  0) 

which  integrates  to  give 

sin  I    (p.-pp)      5(2(t).-2(i)    ) 
(3.26)  <1|     ril2>    =— ^-^  pp/ 

sin  2  (p-]^-p2)       12 

an  expression  which  is  not  normalizable  but  which  shows 
the  typical  strong  anisotropy  of  diametric  correlation. 
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k.      Oscillator  Based  Density  Matrix. 

There  Is  one  Independent-particle  model  In  which 
averaging  out  over  a  somewhat  unusual  broken  symmetry 
has  precisely  the  effect  of  Introducing  pair  Interactions 
Proceeding  In  reverse  order,  one  can  turn  on  harmonic 
Interactions  between  harmonically  bound  particles  by  the 
transformation  of  the  Hamiltonlan 

by  shifting  all  the  coordinates  with  respect  to  the 
center  of  mass 

(4.2)   X.  -^  X.-  A  >   X.  and  letting  K  =  (2> -A  N)ma3  , 


to  give 

2 


2 


Pi  ,  imj'^-K  ^- —  2        K 


1  1  If  J 

In  terms  of  center  of  mass  and  relative  coordinates 
and  momentum 

(4.4)     p^  =  p.  +_^,    P^^IZPi, 

^1  =  ^1  ^  ^  '    ^  =  ^  ^  ^1  • 
The  relative  part  of  Hq  under  our  transformation 
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Is   unchanged,    but   the   center   of  mass   part 
(4.6)  H  =  :^  +  iH^  x^ 

is   transformed   into 

(4   7)  H  -  P       I    11103     -NK      2 

^^''  '  "CM2m+  2  ^• 

Thus  the  corresponding  wave  functions  are  connected  by 
the  unitary  transformation  which  takes 

(4.8)  1  (£?  +"1^  x^) 

^         '  00  ^2m  ^  2  ^    ' 


Into 


2     2 
-       ,jJ_  ,  moo 

(^2_  NK)l/2  ^2m  +    2 


1       (2^^ma4^,2) 


m 

(where  we  have  divided  by  frequency  to  identify  spectra) 
Our  transformation  therefore  takes 

(4.9)  P  into  _^-^P 

moo 

and 

^   .  ^    f^      NK  ,1/4 
X  into   (1 p)    X 

moo 
so  that  a  general  function  F  goes  to  UPU  where 


-  |[l-(l-NK/ma)2)  ^/^]X-P 


(4.10)      U  =  :  e  * 


which  can  be  written  in  the  form 
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'"•^^'"-^u.aMiiA, 


-ft  ^"    ^"  2' 


.  l'K/(l-(l-NK/nn)2)l/^)s.t   is-x  it-p 
e  e  e  ds   dt 


Hence   setting   A  =   l-(l-NK/im)    )    '     ,    we   ha 


ve 


n     Mj^       f^r    i  -^r  s  •  t  is  '  X  .      it  •  p  . 

(2^.12)  U  =  ^l^      //  e       ^  TT(e  Je  J)dsdt 


The  many-body  center  of  mass-magnifying  transformation  U 

is  thus  an  average  of  one-body  unitary  transformations,  any 

of  which  destroys  the  x-p  reciprocity  symmetry  (Born-Yukawa) 

of  the  Hamiltonian.   According  to  our  approximation  method, 

the  basic  generating  transformation  is  given  by 

P    i  ^  s-t  is-x  ^  -is-y 

<x|r-Lly>=/e       ^  ZZ  e  u^(x  +fit)   u^(y+ftt)e  dsdt 

and  r p  is  similarly  obtained.   In  Section  5,  we  carry  out 
a  direct  application  of  the  many-body  transformation  on 
the  N  particle  density  matrix  and  obtain  the  expressions 
for  the  one  and  two  body  density  matrices,  to  within  our 
approximation,  in  terms  of  the  uncorrelated  density  matrices. 
On  the,  basis  of  these  representations  of  p,  and  f~p  some 
properties  of  the  density  matrices  are  easily  shown  and 
application  is  made  of  the  formalism. 

We  thus  have  an  interaction  producing  transformation 
which  we  know  is  effective  for  harmonic  oscillators.   We 
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still  must  obtain  the  basic  Fq-   Consider  closed-shell 
atoms.   To  within  a  simple  coordinate  transformation, 
harmonic  oscillator  and  Coulomb  wave  functions  are  quite 
similar.   If  n   denotes  the  radial  quantum  number,  for 


the  harmonic  oscillator 


1/2 


fu    ih\       I  /o      -^+3/2/  r*         X    /      £  2 

''r^"'  (n^+^+|): 

i+l/2      ^   o        m 
•    L  (ax)    Y.(0,<t))    , 

^r  ^ 

and 

E      „      =    (2n     +  i   +  -2)  -i-T  CD 
n  ^m  r  2 

while  for  the  Coulomb  case, 

U^/2  _   If   ^o   >, 

(4.15)   ^„  z;^  =  P+^/p   ^  ®     ^ 

V""   (n^+^+l)^+^/^ 

(2^+1)    k         m 

•  ^n     ^TTTITl^  )  Y^(«'^) 
r       r 


12  2 
7^  a  X 


and 


V^ 


V"^     4(n^+i+l)2  ' 

but  they  can  be  made  to  correspond  exactly  under  a  one  body 
unitary  transformation.   However  the  Coulomb  level  ordering 
with  energy  is  not  the  same  as  that  for  the  harmonic 
oscillator,  so  that  transforming  the  relatively  easy 
harmonic  oscillator  Fermi  sphere  Pq  ^^  ^'-'^  sufficient. 
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Instead  we  note  that  for  the  Coulomb  Perml  sphere, 

n  4-^+1  <  n.   For  harmonic  oscillators  this  means  that 
r     — 

E/to  <  2n-(^+l/2),  or  L^/h^  <  (2n-  E/Tku)^  -  l/4 . 

The  annihilation  and  creation  operators  for  a  harmonic 

oscillator  are 

(4.16a)   a  =  ^ p  +  1  /^  x  ,   b  =  a*=  -J^x-1 /^ 

or 


(4.l6b)   p  =  /^  (a+b)  ,    X  =  / 


■f\       a-b 


2ma)   1 


and  then 

1    ^ 
(4.l6c)   ^  =  — ^  b  ^Q 

"    yn  I 

for  each  dimension. 
Now 

(4.17a)  L^  =  x^Pg  -  X2P3_   =  y(^i^2"^2^1^ 

so 

(4.17b)  L  =  -h/l  b  X  a 

and 

(4.17c)  E  =   (b  •  a  +  |>tn)  . 


2 


Further, 


(4.18)   (b  X  a)^  +  (b-a)^  =  : (bxa)^:  -2b-a  +  :(b.a)^:+b-a 

2  2 

=  b  a  -  b  •  a 
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so  that  the  closed  shell  condition  becomes 

(4.19)  0  <  (.fl  -  |)%  (2n-  |)^  -  2(2n-  |)  (jf^  -ij-^-^ 


=  b^a^  -  (^  -  ^)    2(2n-l)  +  4n^  -  6n  +  2  . 

This  yields  a  one  body  density  matrix  (with  Z  a  unit  step 
up  at  the  origin)  of 

(4.20)  \~Q  =  -RZih^a^-   2(2n-l)(b  •  a+ |)  +  (2n+l)  (2n-l)  ]  , 

which  of  course  can  be  used  for  annihilation-creation 
structures  other  than  those  of  the  oscillator. 

When  we  apply  the  P -,  approximation  method  to  a 
pair  of  oscillators,  in  Section  5,  the  result  is  found  to 
be  exact.   This  is  not  far  from  true  for  any  set  of 
Independent  oscillator  states.   The  reason  is  that  due  to 
Eq.  (4.1 6b)  the  unitary  transformations  required  are 
very  special.   They  are  exponentials  of  linear  combinations 
of  annlhllators  and  creators 

(7-a*-7*-a) 

(4.21)  U   =  e 

7 

with  the  property  of  having  a  closed  form  inverse  on  any 
N  by  N  subspace  without  holes.   By  a  subspace  of  orbltals 
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without  holes,  we  mean  one  taken  into  itself  by  a. 
It  then  follows  that  for  the  states  p,q  <  N, 


*  * 


pa     aa  oo  pa  aa 

(^.22)     <p|e        e       |  q>     =  :^L_  <P  I  e         |^i><M-|e       |q> 

N  pa*  aa 

=  ZZ<Cp|e         lM-><;'^i|e      |q^ 

while 

N  -aa  aa  oo  -aa  aa 

(4.25a)     ZZ<Cp|e         ll^><l^|e      |q>  =  ZZ<Ple         lM-><U|e      |q^ 
M.=  l  M.=  l 

pq 


with  its    conjugate 

N  p 

(4.23TD)      ^<q|e         |p.><^i|e'         |p>     =5 


N  pa*  -pa* 


H=l  P^ 


Hence   if 


(4.24a)  U        =  <u    I    U    lu    > 

pq        ^   p  '  '    q^ 


we   have 

1  *  -x-  * 

-    p   7-7  7B.        -7   a 

(4.24b)      U        =    e      ^  <p|e  e  |q> 


or 


1         *  *  * 

--2  7-7      N  73*  -7   a 

;  XZ'sPle         lM-><l^|e  |q> 


_-Y  2  y -y        N  7  a  -7-a 

(U      )         =    e  XZ<^|e         IM-Xtile  |pS 

^P  M^l 


while 
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N   _.    .  *  * 


(4.25)      Det   U  =    e      ^  Det<q|e  |tx;>    Det  <tx  |  e        ^|p> 


{y-7'){i-y')* 


N  ^ 

=   e 


Thus 

_    N^ 
(4.26)      Det<u    |U^,U    |u   >     =    e        ^ 

UT^P  P        /        /       u 

•  e  ZI<^I^  I^L><^L|e  |p>. 

|X=1  '^ 

In   our   approximation  we  require   the  replacement   of 

by 

-  |(77^+7'7'*-27'*-7)  ^      7a''-7'''a  -7a*+7'*a 

^ToU/e  foe  )U^,rc 

with  a  similar  direct  modification  for  pp. 
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5-    Correlation  through  Transformation. 

The  many  body  transformation  of  Section  k   will  here 
be  used  directly  to  Introduce  correlation  and  obtain  the 
approximate  oscillator  based  density  matrix.   The 
coordinate  dependence  of  the  diagonal  elements  will  be 
shown  to  be  exact.   For  the  case  of  two  particles  In  a 
harmonic  oscillator  potential  our  density  matrix  will 
give  the  exact  result  for  the  harmonic  oscillator 
Interaction.   In  the  nonlnteractlng  N  particle  density 
matrix,  formed  by  Slater  determinants,  displace  each 
coordinate  by  a  fraction  of  the  center  of  mass  of  the 
system.   Let 

(5.i)x    =x+x.  ,    x  =  XZ  ^n-  >    y    =  y  +  Yn  >  y=  ZZ  y^  • 

■^        1=2   ^  1=2   ^ 

Our  determinants  will  be  specified  symbolically  by  giving 
one  of  their  elements. 

(5.2)  <...x^...|r^^M-..yi...>  = 

,      _  t       *      _  t 

=  n  Det  u  (x.  +  Ax  )  Det  u  (y.+  Ay  ). 

Within  a  normalization  factor  n  which  we  omit  in  the 
following, 

(5.3)  <x|r^'^^i)|y>=  Tr2^^^^<x|r^^My> 

-ls(R-x) 


dXg  .  .  .  dXj^  /  /  e 
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•  Det  Up(x^+  AR  +  Ax-^)  Det  Up(x  +  XR   +Ay-j_ )  dR  ds; 
expanding  each  determinant  by  minors  this  Is  (letting 

^1  ^  ^'   yi  =  y): 

rr  -is.R  ^ 

=  JJ    dR  ds  e      2^  u^(x  +  AR  +  Ax)  Ug(y  +  AR  +  Ay)  • 

r ,  s 

where  Q_^  =  cofactor  (a   ). 
rs  rs 

,     p    Is-x.   ^ 

In  these  a^p,  let  x^  =  x^+  AR,  so  a^p  =  e"^^^ "^  a    ,  with 

Is  "X, 


(5-5)      ttqp  =  J   Up(x^+  Ay)    e  ^  u^Cx^  +  Ax)    dx^    , 

and  Q        =   cofactor    (a^^ ) . 

,_rT\Ti  rr  -1s-r[i+(n-i)a] 

(5.6)    <x|r^^^(i)|y>  =  j[/  ds  dR  e 

.    •   YZ  u^(x[l+A]+AR)    u*(y[l+A]+AR)    Q^g[s,y,x]. 
r ,  s 

This  expression  is  exact  and  derived  from  a  known  N  particle 
wave  function  with  we  have  reduced  to  a  double  (3  dimensional) 
integral  of  products  of  single  Integrals.   We  can  simplify 
this  further  (--  with  respect  to  getting  a  tractable 
expression)  by  the  following  procedure.   We  note  that 
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Q   =  a   det  a,  but  evaluation  of  a    is  a  laborious 
rs    sr 

undertaking  at  best.   We  may  write 

-  ^  Ayp  Is-x.  ^  Axp 
^5-7)     -qp=<-pl-       e    ^e^    \u^y  . 


When  higher  order  commutation  brackets  give  zero,  we  have 
the  operator  equation 


a(A+B)    aA  aB   -3  a^[A,B] 
(5.8)     e       =  e   e   e 


We  may  therefore  express  the  operator  a  as 

-fAyp   is-x^     ^Axp 
(5-9)     a  =   e  e  e 

-1  ^(x+y)       (-  |Ayp+lsx^+  -^Axp) 
=   e  e 

So  that  for  a  complete  set  (00  matrix)   the  Inverse 
element  Is  Immediately  known 

(5.10)     a^p  =  /^p(^l+  ^^^®      "^  "q^^l^  ^^^  ^^1  • 
00 

We  will  also  use  the  relationship 

y    (In  A) 
Tr  In  A    V       ^^ 
{5.11)  Det  A  =  e        =  e 

We  will  now  use  these  relationships  to  approximate  Q 
by  the  inverse  of  the  00  matrix  multiplying  the  N  by  N 
determinant . 
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-  -^   Axp   -Is-x.  ^Ayp 
(5.12)   Q,,  =<u  |e  ^     e     ^  e^    |u> 


rs   ^  r '  '  s 

> 

1=1 


exp  {-1   I  Ns(x+y)+2Z  ("  :^  Ayp+isx^+ :i  Axp) 


For  states  of  definite  parity  the  sum  which  has  linear  arguments 
in  X  and  h/hx   vanishes.   To  the  above  approximation  we  .,ow  have: 

,       /^J^  ,  rr     -ls.R[l+(N-l)A]       -±2^U+y) 

{5.13)    <x|r^^^(l)ly>=    //   e  e  2 


N,N  ^  n  :^ 

^  u    (x[l+A]+AR)u    (y[l+A]+AR)- /   dx-,    u    (x-,+Ax) 


N,N 
r,s 
-is  -x 


u    (x-,+Ay)    ds   dR 


Taking    /   ds  and   then    /   dx^ , 


'   dR  21  u    (x[l+A]+AR)u    (-R[1+(N-1)A]-  -^    (x+y)    +  Ax) 
Ug(-R[l+(N-1)A]    -  ^    (x+y)    +  Ay)    u* (y[ 1+A]+AR) 


(5.14)        =    TdR    r  (x+A[x-ffi],-R+A[x-m]-  ^(x+y+2R)) 


•r(-R+A(y-H^)    -~   (x+y+2R),    y+A[y4R])    . 

For  A  =  0,    we   obtain 

(5.15)  r(x,y)  =<x|r^^^i)ly> 

as  required   for  a  legitimate  density  matrix  of  noninteracting 
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particles  (i.e.  T  f (x,R)  f (R,y)  dR  =  r(x,y)).   The 
second  order  density  matrix  Is  (without  the  normalization 
factor  n) 

(5.16)  <X3_X2I  r  ^^^  (2)  Iyiy2>=  Tr^  _  _  _^<.  .  .x.  .  .  .  |  f  ^^^  l-yi-) 

A.  orN   -lS-(R-Xp) 

=  /  dx   . . .  dx^   //  e  dR  ds  • 

•  Det  u  (x^+AR+Alx^+x^))  Det  u  (x  +AR+A(y^+y2) )  , 

N 
with  Xp  =  >   X..   We  now  use  a  "double  minor  expansion" 

1=3 
on  each  determinant  to  separate  nonlntegrated  terms. 

Analogously  to  Q 

(5.1Ta)   Q^3  =  cofactor  (a^ga^^)  =  (ag^a;;;^  -  cc'^^^^)    c3et  a  , 
tu 

which  vanishes  for  r  =  t  or  s  =  u,  with 

r    *  Is-x. 

(5.17b)   a^p  =  /  Up(x^+A[y^+y2])e      Uq(x^+A[x^+X2] )  dx^  . 

Now  (5.16)  is 

PP  -ls-R[l+(N-2)A] 

(5.18)         //  ds  dR  e 

N  * 

•  ZZ   u^(x^[l+A]+A[R+X2])Ug(y^[l+A]+A[R+y2])- 
r, s, t,u 


u^(x2[l+A]+A[R+x^])u^(y2[l+A]+A[R+y^] 

Q^s(s,y-^+y2,    x^+X2)     • 
tu 
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We  now  follow  the  same  procedure  as  for  PCD  so  for 

states  of  definite  parity  we  obtain 

.  AN    /„, 
-1  -p-  s-(X)  p    :^  _      -is-x. 

(5-19)   Qrs  =  e  [  /  u^(x^+Ax)e     "^  u^  (x^+Ay)dx^' 

tu 


^       _   -is-x. 

u,  (x.  +  Ax)  e     "^  u  (x  .  +  Ay)  dx  .  - 
t'  J  u'  J     -^     J 


/*       _    ~is  •  X .  _ 

u  (x.  +  Ax)  e       u  (x.  +  Ay)  dx.  • 
r   1  u'  1    '^     1 

/w  —is ■ X .  _ 

u.  (x.  +  Ax)  e     -^  u^(x.  +  Ay)  dx  ]  , 

with  X  =  x-j^+x^,  y  =  y]_+y2>  ^  =  ^+yi  cCq-^   has  been  used  here. 

oo 
Put  this  in  equation  (5.I8),  and  we  obtain  5 (x.+x .-ffi(l+(N-2) A]+ 

+  (AN/2)X).   Let  ^(A,N)  =^=   1  +  (N-2)A.   We  obtain 

(5.20)  <^x^x^\r^^\2)\Y^y^'^=  J  6rJ  6l<ri^:^+HR+x],l+^x)    . 

•  r(x2+A[R4-x],-x  -RX-  ^  X+A5^)-[r(x+Ay,y^+A[R+y]). 
"  •  r(-x  -R^-  ^  X  +  Ay,y2+A[R+y])-r(x+Ay,y2+A[R+y])- 

.  r  (-X  -  R^-  ^  X  +  Ay,  y^+  A[R+y])]>  . 

The  correctness  of  the  functional  form  obtained  for 

diagonal  matrix  elements  by  our  formalism  is  proven  as 
follows: 

For  X  =  y,  let  x.  =  x.  +  Ax,  so  that 
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t 
-IS'AX      n      Is-x.         ^         ,  , 

(5.21)      a        =    e  /    e  '^   u^(x.)    u    (x.  )    dx. 

^■^  qp  J  piqi  1 

-is -Ax 

Cqp(s)     . 

Let  Q,      (s)    =    cofactor    (c      );    so   that  now 
rs  rs 


(5.22)  <x|  r  ^^^  (1)  |x>  -    /7  dR  ds   e 


■ls-R[l+(N-l)A] 


-is-x   A(N-l)  ^ 

•    e  y~  u^Jx[l+A]+AR)u„(x[l+A]+7vR)Q^„(c); 

X  o  X  o 

r ,  s 

_1  Is -Ax       n      ^  -ISX, 

(5.23)  a        =   e  /   u    (xt  )    e  u    (x-,  )    dx-, 
\ji     ^1      v^qp                        J      p'    1'  q'    1  1 

GO  .  --  -     - 

lAs-x 
=    e  aqp(s)    ; 

_1  -iANs'X 

(5.24)  Q^^(x,x)    =  ctqp(x,x)    e 

00  00 

-lAs-x(N-l) 

aqp(s)     . 

But 

-is-xA(N-l)    - 

(5.25)  Q^3(x  x)    =    e  ^s^^^     • 

So  our  approximation  of  the  inverse  by  the  operator 

inverse  in  a  complete  basis,  in  the  case  of  diagonal 

terms,  gives  the  exact  x  dependence  and  only  approximates 

the  s  dependence. 

To  apply  our  formalism  to  the  harmonic  oscillator  in 

the  case  of  two  particles,  the  initial  density  matrix  is 

-an(x2+y2) 


(5.26)       r(x|y)  =  e 


^0' 


36 


Using  this  in  Eq .  (5.14),  we  obtain 

(5-27)    e  °  e  °   2 

P     -a  [R^-2(A^+c^)+2AR-(x+y)(L+c[N-l])] 
•  /  dR  e   ° 


with  L  =  1+A,  c  =  1+(N-1)A.  For  N  =  2,  c  =  L. 
Note  that  at  N  =  2,  the  x-y  cross  term  which  would  make 
this  expression  a  non  legitimate  density  matrix  vanishes. 
For  N  =  2  carrying  out  the  R  integration  gives 

(5.28)  e  °  ( 1^-^)5/2  e  o  a2+L^ 

2aQ(A^+L^) 

which  is  a  convergent  quantity  for  all  values  of  A. 
Let 

(5.29)  A  =  aQ(L^  +  A^)   and  B  =  Qq  4aL  . 
For  our  case  here 

(5. 30)  r^^^  =  /(x^  Xg)-  ^(y^  y^) 

-A(x2+x|)-Bx-,_-X2   -A(y^+y2)-By^-y2 


by  inspection  of  (5.28)  when  N  =  2.   This  however  is  exactly 
the  harmonic  oscillator  interacting  wave  function.   This 
result  which  yields  the  exact  form  of  the  density  matrix  for 
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the  interacting  harmonic  oscillator  case,  to  within  the 
assignment  of  the  parameters,  is  very  unusual. 

Aside  from  our  approximation,  in  deriving  our  final 
formula,  of  the  matrix  element  which  should  be  good  for 
large  systems  but  poor  for  small  systems   the  entire 
derivation  assumed  a  determinantial  form.   For  the  N  =  2 
ground  state  due  to  separation  of  space  and  spin  parts, 
the  separate  space  part  which  is  symmetric  certainly  does 
not  satisfy  this  derivation.   It  therefore  seems  that  the 
final  equation  has  a  wider  validity  than  its  derivation 
indicates.   That  this  is  so  has  been  formally  shown  in  the 
last  part  of  Section  k. 

We  now  wish  to  deal  with  atoms.   The  harmonic  oscillator 
solutions,  as  they  stand,  have  quadratic  arguments  and  so 

do  not  give  reasonable  approximations  to  atomic  states. 

2 
However  we  may  carry  out  the  operation  of  letting  the  r^ 

be  transformed  to  linear  arguments  p..   This  transformation 

carries  the  noninteracting  harmonic  oscillator  solution 

for  two  particles  into  the  two  particle  noninteracting 

-bCpn+pp) 
coulomb  solution  f  =   e  .We  note  also  that 

Schwinger  has  thus  reduced  the  radial  equation  of  the 

Kepler  problem  to  the  radial  equation  of  the  harmonic 

oscillator  [10].   This  leads  us  to  expect  that  the  same 

transformation  applied  in  conjunction  with  our  formalism 

will  yield  good  correlation  energy  corrections  for  atomic 
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systems.   We  carry  this  out  by  applying  Eq .  (5.l4)  to  the 
harmonic  oscillator  case,  and  then  carry  out  the 

; A 

r.  =  yp .  p.  transformation.   Note  that  this  Is  the  case 
where  our  formalism  gives  exact  results  so  our  result  Is 
Identical  with  that  which  would  be  obtained  by  operating 
directly  on  the  interacting  harmonic  oscillator  solution. 
The  resulting  trial  wave  function  is  therefore  determined 
to  be  (by  Eq .  (5-30) ) 

r.  -.1^        ,    M   -  S(Pl+P2^-  I/P1P2  Pl-P2 
(5. 31)       ^  =  N  e 

Details  of  the  calculation  are  presented  in  Appendix  I, 
The  value  we  obtain  for  the  ground  state  energy  is 
E  =  -2.86765  (for  Qq  =  3-2520  and  A  =  0.02014),  which  Is 
very  close  to  the  value  obtained  by  the  "harmonic  atom" 
method  in  Section  6,    and  shows  that  our  approach  accounts 
for  the  correlation  energy  very  well . 

The  direct  application  of  our  formalism  to 
hydrogen-like  orbitals  for  N  =  2,  using  the  ground  state 
non-interacting  one-body  density  matrix, 

(5-32)  r  =  e  " 

yields 


(5-33)  <x|  I  (l)|y)>=  /  dR  exp  ^ -a  (  \/x^L^+  A^R^+2x-RAL 

+  /'R^c^+cAR-(-2x+N(x+y))+A^x^+(^)^(x+y)2_NA^x-rx+y) 
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•  exp  <(-  q(  >/y^L^+AS^+2y-RAL 


+  yR^c^+cAR.(-2y+N(x+y))+AV+(^)^(x+y)^-NA2y(x+y)y; 


set  N  =  2, 


(5.3^)   =  /  dR  exp  < -QqC  yx^L^+A^R^+2x-RAL 


/r  L   +  2ALR-y  +  Ay)  - 


exp  \-cLq{   >/y^L^+A^R^+2ykAL  +  yR^L^+2Am  -x+A^x^) 


dR  exp  ^-2aQ(/xL^+A^R^+2ALx-R  +  yy^L^+A^R^+2ALy -R  ) 


We  note  in  this  case  also  that  N  =  2  has  the  special  property 
—  needed  for  a  legitimate  N  =  2  density  matrix,  but  not 
needed  for  N  >  2,  that  the  x-y  terms  cancel  out. 

Again,  by  inspection  of  (5.3^)  if  we  let  R  ^^ Xp,yp 
we  can  form  a  f    (x,Xo)  ^(y-iYp)  which  will  yield  this 
interacting  density  matrix.   The  wave  function  is 


(5-35)   ^(x^Xp)  =   exp  <-   OqC  /x^ 


2   2  2     -   - 
L  +A  X2+2ALx^-X2  + 


+  /x2L^+A^x^+2ALx^-X2 


(neglecting  normalization),  which  also  satisfies  the 
symmetry  requirements. 
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This  form  of  f^due  to  the  angular  argument  under  the 
square  root  sign, gives  rise  to  difficult  Integrals,  and 
we  shall  not  use  It  for  computations. 

The  transformation  we  have  introduced  gives  an  exact 
solution  for  the  interacting  harmonic  oscillator  case. 
It  is  reasonable  to  suppose  that  in  the  general  case  a  sum 
of  such  solutions  may  give  a  reasonable  wave  function,  l.e, 

(5.36)   /  -  ^  g*(\)  Det  Up(y^+  Aj^y'). 

¥e  may  consider  A  a  continuous  parameter  so  that 

(5-37)    ^(x^...Xjj)  =  rg(A)  Det  Up(x^+  Ax')  dA  . 

This  is  recognized  as  a  form  of  Eq.  (3.7).   Here  we  shall 
proceed  by  direct  means  to  obtain  the  one  and  two  particle 
density  matrices  in  terms  of  the  initial  nonlnteracting 
density  matrices. 

Let  us  denote  the  N-th  order  density  matrix  of  the 
system  by 

(5.38)  <..-x^...|  r!^M...yi>   =  j^J'dA  dA'  g(A)  g*(A'). 
•  Det  u^(x^+Ax')  Det  Up(y^4-A'y') 

=  TdA  dA'  g(A)  g^A')  <...x.  ...|r^^^|...y.  •••> 
J  ^      AA'      1 

and  the  first  order  density  matrix  Is  (to  within  a 
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normalization  factor) 


=  J^dA  dA'  g(A)  g'^CA')  Tr2___j^<.  ..x^|T^^||  ...y^...> 

=  Jd-x  dA'  g(A)  g*(v)<x|r^^!(i)|yX 

so  that  our  previous  formalism  can  be  carried  out  and 
at  the  end  we  apply  /  dA  dA '  g(A)  g  (A')  to  the  result. 
However  we  must  make  some  slight  modifications  due  to  the 
presence  of  A  and  A'. 

(5.40)<x|  r^^^(l)|y>=  jJdA  dA'  g(A)  g*(A')- 

pp  -is  "R  ^ 

•  //  dR  ds  e      XI  u^(x+AR+Ax)Ug(y+A'R+A'y)  • 


Qjlg[s,A'(R+y),A(R+x)]  ,   with  Q^g  =  cofactor  (a^^) 


(5-^1)   a'   =  /"up(x^+A'[R+y])u  (x^+A[R+x])  dx^  e 


(due  to  A',  R  remains  in  Q'); 

-^A'(y-m)p   is-x.   i  A(x-f^)p 
-<Up|e  ^         e    ^  e^        |Uq>. 

So  the  operator 
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-  ^(A'(y-ffi)+A{x+R) 
(5.^2)   a   =   e       "^ 

1^. 


{-  -pjA'(y4^)p+isx^+ ^A(x-m)p} 


For      states      of       definite  parity, 

_.•   _Ns^^  l„■^„^      INs 
(5.^3)      Det  a   =   e 


-1-f  (A'y+Ax)-^R(A+A') 


-1  p      =f,  -IS'X. 

(5.^^)      a^p   =  j   ^p(Xi+  A{R+x)    e  ^  Uq(xi+  ^'(y-f^))    dx^ 


00 

Now  we  use 


(5.^5)  Q^3  =  a-J:  det  a 

00 

N, 


LetA=      1   +|(A+A'),   B-  |(^x+A'y).      Then 

(5.46)    <x|r!^^(i)|y>=  j^TdA  dA'   g(A)   g''(A')- 

po  -lS'(R,4  +  S)    M  ^ 

•    //   dR   ds   e  XI  u    (x[H-A]+AR)u    (y[  1+A  '  ] +A 'R) 

^^  r,s      ^  ^ 

-is -Xn 


Take 

and   get 


•  J   dx^  u^(x^+A[R+x])    e  ^  "s^^l"*"  ^'t^+y]) 

/  ds  -^     5(Ar   +\^+x^).      Then  take    /   dx. 


(5.^7)     j/JdA   dA'    g(A)    g(A')     fdR  T~  u^(x[l+Al+AR)- 

r ,  s 

•    u*(-[/i^R+B]+A[R+x])u^(-[AR+6]+A'[R+y])u^(y[l+A']+A'R)  = 
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=  /y  dA  dA'  g(A)  g(A')  /  dR  r(x+A[x-m],-[AR+B] 

+  A[R+x])  r(-[AR+B]+A'[R+y],y+A'[y-m]). 

The  expression  for  the  two  body  density  matrix  can  be 
obtained  In  exactly  the  same  manner. 
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6.   The  Harmonic  Atom. 

Since  the  Intrinsic  state  approach  can  be  used 
without  approximation  for  harmonic  Internal  and  external 
potentials  by  working  with  orbital  annlhllators  and 
creators,  the  same  would  hold  for  any  set  of  orbltals. 
This  Is  Identical  to  choosing  as  an  ansatz  a  one  body 
unitary  transform  of  a  totally  harmonic  atom.   Such  an 
ansatz  clearly  has  an  Invariant  characterization 
consisting  of  the  eigenvalues  of  the  one  body  density 
matrix  and  the  full  two  body  density  matrix  In  the  T  -, 
representation.   To  see  if  these  are  reasonable  for 
application  to  atoms   let  us  first  look  at  the  two  body 
problem.   For  a  two-electron  system,  with  opposite  spin, 
the  space  wave  function  is  symmetric  and  we  may  define 
the  symmetric  matrix 

(6.1)  <l|^|2>=  ^(1,2)  . 

Hence  the   one  body  density  matrix  is  given  by 

(6.2)  r]_  =  2^^*. 

If  also  Tp   is  real  then  the  real  orthogonal  transformation 
which  diagonalizes  f   diagonalizes  \~ -t  >    and  this  is  our  case. 
Now 

IT 


^5 


for  the  interacting  ground  state.   J.  K.  Percus  has  shown 

that  the  eigenvalues  of  the  density  matrix  go  as    ■   " 

(n^+n2+n^)c 
c^e        -^       and  I  have  checked  that  the  results  obtained 

are  quite  close  to  those  for  the  helium  calculation  by 

Auffray  and  SobolfS].  \^ ^=-^{\,2)   ^*(1'2')  In  the  same 

representation  does  not  have  the  correct  sign  for  all 

elements,  but  of  course  does  for  the  crucial  diagonal  ones. 

The  method  thus  should  give  reasonable  results.   A  possible 

approach  to  the  N  body  problem  with  Interaction  thus  Is 

available  to  us.   The  nonint era c ting  density  matrix  for  N 

particles  in  a  harmonic  well  as  given  by  Eq.  (4.20)  can 

by  the  formalism  in  Eq .  (4.27)  be  transformed  to  the 

interacting  harmonic  oscillator  density  matrix.   We  would 

then  desire  to  use  unitary  transformations  to  a  more 

appropriate  set  of  orbitals.  We  shall  now  study  this 

problem  in  a  more  direct  manner,  and  its  application. 

Consider  a  system  for  which  the  solution  ())(1...N) 

is  associated  with  the  solution  for  another  physical  system 

^(l'...N')  by  means  of  the  unitary  transformation 

U(l. . .N,l' . . .N') . 

(6.4)   (t)(l  ...N)  =  Tud.  ..N,l'.  ..N')^(l'.  ..N')  dXj^. 

The  density  matrix  for  (f)  will  then  be  given  by 
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(6.5)  r(|,(l'2'|l  2)  =  JJj   /(l"...N")U*(l"...N"|l',2',3'"-.) 

U(l,2,3'"  .  .  .N"'  |l",2^  .  .N"')^(1'^  .  .N'")  (dx")  (dx'")  (dx/J^')  . 

This  would  represent  an  exact  solution.   On  the  other 

hand,  reducing  our  U's  to  product  U's,  and  our  tj/^s    to 

antlsynunetrlzed  products,  we  would  obtain  a  solution  of 

the  Slater  determinant  type  In  which  electronic  correlation 

Is  neglected.   We  will  Introduce  correlation  by  retaining 

the  model  wave  function  in  its  Interacting  form  but  assuming 

that  a  reasonable  approximation  to  the  true  solution  can 

be  obtained  using  U's  of  single  particle  forms  i.e. 

N 
{6.6)  U(1...n|1'.  .  .N')  =   ^T  U(i,i')  . 

i,i'=l 

With  this  assumption, 

(6.7)  y'u*(l°...N°|l^,2^,3'.-.N')U(l'...N' |1  2  5...N)(dx^,2,) 

-  U*(1°1^)U*(2°2-^)U(1'1)U(2'2)  Tl      5(1^-1)  • 

i/1,2 

So 

(6.8)  ri(l'2'|l  2)  =  ru''(l°l')U*(2°2') 

•  r^(l°2°|l"2")U(l  1")U(2  2")  dl°  d2°  dl"  d2" 
and 

(6.9)  E,  =   r  r^,(l°2°|l"2")[ru''(l°l')U*(2°2') 

•H(l  2)U(1  1")U(2  2")  dl  d2]  dl°  d2°  dl"  d2"  , 
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where 


(6.10)     H(l  2)  =  N  H(l)  +  N(N-l)  ^(3_  2) 


H(l)  are  the  one  body  terms  of  the  Hamlltonlan  and  V(l  2) 
are  the  Interaction  terms. 

The  simplest  application  of  our  formalism  Is  to  let 
our  unitary  transformation  be 

(6.11)  U(l,l')  =  6(1-1')  . 

Applying  this  to  the  wave  function  of  the  central  harmonic 
oscillator  potential  with  harmonic  oscillator  interactions 
between  the  two  particles  will  yield  a  poor  estimate  of  the 
energy  but  will  give  a  clear  indication  of  the  dependence 
of  the  energy  on  the  Interparticle  Interaction  of  our  model. 
The  wave  function  used  is 

-  — (r   +r    )+Br    -r  2      2/ 

(6.12)  ^  =  N^,.  -  N  e      2      ^     ^  ^      ^     where  N^=    (A_i:|_)5/2 


as 


IT 


~    ^    -A(r2+r2) 
e 


(6.15)  N^  r^y^y  di  d2  =  N^(47r)^  fr 


0   0 

slnh  2Br, rp 

r-,r^  dr,    dr^ 


2B  "1"2        1        2 

/       S°  -(A-BVA)r^ 

=  N^   47r(|)^/2  j   r^  e  ^   dr^     =      1    . 

0 
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(6.14)   E  =  l^f*{l   2)H(1  2)^(1  2)  dl  d2 


H(l  2)  =  -  ^  (V?+v2)  _  2|!  .  2|!  ^  e2 


2m  '  1   2     r^    r^        "^i   2 


Considering  the  kinetic  energy,  potential  energy  and 
interaction  energy  separately  we  obtain 

(6.15)     g    /,^(.v2.v2),  ,id2=^.f      . 


r/(-^  -^)^  dl   d2  =   -   kZe^    / 
J                   ^1  ^2 

r  ,*^    e^^,  ,.     ,o  2    /2(A+B) 

/    ^    (t^ — )^  dl   d2  =   e    V      ^         ' 

J  ^2.2  ^ 


ttA 


S  o   in  a  .  u  . 

(6,16)       E  =  f  -  8  yAp?    ^  y-EMl  . 

In  terms  of  the  strength  of  the  central  potential  k 
and  the  interacting  potential  K, 


(6.17)  E 


=  i  [  v/kT2K  +    v/k   ]    -   4Z    /I^^Kl^     + 
^  V  7r(v/k+2K+\/k   ) 


^        IT 

For  K  =  0,  setting  ^E/Sk  =  0,  we  find  the  minimum  for 
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(6.18)  vS:  =  I  (  -^ ^^  )   where  a  =  \/k  , 

^   \/27r     v/ir 

so  E  =  -  2.3002  a.u.  for  k  ^   2.353.  •  •- 

If  K  is  positive  and  increases  the  energy  gets  larger. 

If  K  is  negative  and  |k|  increases  the  energy  decreases. 

The  minimum  of  -2.32389  a.u.  occurs  at 
k  =  2.8680,  and  K  =  -  0.4460. 

A  plot  of  the  variation  of  the  energy  versus  K  with 
k  fixed  yields  a  parabolic  form  and  the  same  is  true  for 
the  energy  versus  B  with  A  fixed. 

We  see  that  we  get  an  improvement  in  the  energy  of 
-O.O23  a.u.   The  main  point  is  that  the  variation  in  energy 
follows  the  model  in  the  sense  that  when  K  is  positive, 
representing  an  attractive  interaction  between  the  particles 
which  does  not  correspond  to  the  helium  atom, the  correlation 
introduced  increases  the  energy.   However  when  K  becomes 
negative  the  energy  is  lowered  by  the  correlation 
Introduced  as  it  corresponds  more  closely  to  the  correct 
correlation  effect  for  Helium.   The  overall  energy  value 
we  have  obtained  is  very  poor  due  to  our  using  the  harmonic 
oscillator  wave  function  directly.   We  therefore  must  use 
a  more  reasonable  unitary  transformation  U(l,l'). 

When  the  simple  product  model 

-b(r  +r  ) 

(6.19)  ^H^p_  =  N  e    ^   2   ^   Vq{1)vq(2)  , 
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which  is  a  solution  for  the  helium  atom  without  electronic 
interaction  is  used  to  obtain  the  energy,  the  minimum 
value  obtained  is  E  =  -2.847656O  a.u.  for  b  =  27/16  = 
1.68750,  which  is  a  very  good  result  as  the  true  energy 
is  -2.9037  a.u.   We  therefore  seek  a  U(l,l')  such  that 
the  correlated  wave  function  it  produces  will  reduce  to 
tpy,  -p     as  the  electronic  interaction  goes  to  zero  in  our 
model  i.e.  ■^tj  q  -^  Uq(1)uq(2)   a  product  form. 
Ideally  we  would  use 

00         ^ 

(6.20)  U(l'l)  =  XZ   v,(l')u,(l)  . 

i=0   ^     ^ 

As  this  form  is  not  known  to  us  we  will  use  a  U  which  will 
meet  our  specifications  with  respect  to  the  ground  state 
but  will  connect  modified  sets  for  higher  states.   Such  a 
U  is 

u  (l)-v^(l)      ^      ^ 

(6.21)  U(l,l')  =   (1-1-) ^ ^- [uq(1')-Vq(1')] 

l-|^UQ(2)vQ(2)d2 

A  detailed  calculation  using  this  U(l  1')  follows  in  the 
application  of  the  formalism  to  the  calculation  of  the 
ground  state  energy  of  helium. 

The  integral  expressions  for  the  energy,  E,  can  be 
set  up  in  two  equivalent  forms  in  our  formalism. 
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(6.22)  E  =  [  i}j*{l2)[     ru''(l"l)U*(2"2)  H(l  2)U(l  l')U(2  2')dld2] 


•  ^(1'2')  dl '  d2"  d2'  dl" 

=    N(N-l)  /  r^(l"2"|l'2')[y^U^(l"l)U*(2"2)H(l  2). 

•  U(l  1')  U(2  2')  dl  d2]  dl'  d2 '  dl"  d2"  . 

(6.25)   E  =  f  ^*  [1   2)  H(l  2)  (t)(l  2)  dl  d2 

with 

Ml  2)  =  Tud  1')  U(2  2')  ?//(l'  2')  dl'  d2'. 

The  single  particle  product  wave  function  to  be  used  as 
our  approximation  for  the  helium  wave  function  in  construct- 
ing our  single  particle  unitary  operators  U(l  1')  is 

3  -b(r.+rp) 

The  noninteracting  hannonic  oscillator  wave  function  for  two 

particles 

_  a/  2   2, 
(6.25)     -P^^Q^    =   Uq(1)  Uq(2)  =  (^)5/2  J   2  ^1  ^2 

is  also  used  to  construct  U  in  such  a  manner  that  if  ^^^  ^ 
is  substituted  for  ^  in  (6.2  2)  (|)  in  (6.2  3)  becomes  ^„  „  . 
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To  Improve  our  energy  estimate  we  will  operate  on  an 
"Interacting"  wave  function  to  take  care  of  the  effect  of 
correlation  between  the  electrons.   We  can  denote  the 
unitary  form  (6.21)  by 

(6.26)  U(l  1')  =    5(1-1')  -  A(l  1') 

-  5(1-1')  -_R(1)  ^  JlTd')    . 


Note  in  our  case  the  quantities  are  real,  so  u„  =  u^, 
^0  =  ^0- 


n  u  (i)-vq(i) 

(6.27)  Ji-d)  °  ° 


/       } 

V   1    -    /    Uq    Vq    d2 


where 


00       a   2  , 


(6.28)  Juo(l)Vo(l)    dl=    (^)5/\^)l/2.w|e'    2    1'      1   r^d 

0 

5/4,  3/2  ,  I -b^/2a    H         ^2  , 

TT  '  a  ^  a  v2a 

The  interacting  wave  function  to  be  used  is  a  solution  of 

the  interacting  harmonic  oscillator  problem  with  the 

Hamiltonian 

2 

(6.29)  H^'ZZ  (ii+|x^)  +|ZI(xi--)'  . 


where  we  shall  also  have  K  take  on  negative  values. 


^1 
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The  symmetric  ground  state  solution  Is 


(6.30)     ^j^  =  N  e 


a 


Xj^+  2j;j(p-a)(ll  x^; 


with  a   =  miu/fi  =    N/ink/fi      (k=iiuo^),    p   =  \/m(k+NK)/fi 
For  N  =   2, 

^^2_g2   3A      -  I  A(r^+r^)+B(?^-?2) 


(6.31)  ^^  =    (- 


TT 


with  A  =  -|(a  +  p),    B  =  ^(p  -a) 


The  wave  function,  ij/y-,    Is  now  used  with  H(l  2)  for 
helium  In  Eq.  (6.23).   The  details  of  the  calculation  are 
given  In  Appendix  II.     Part  of  the  Integration  Involved 
was  carried  out  numerically  on  the  computer.   Varying  the 
parameters  led  to  a  minimum  energy  of  -2.868  a.u.  for 
k  =  1.4  and  K  =  -O.13,  representing  an  Improvement  of 
0.02  a.u.  over  the  uncorrelated  wave  function  with  K  =  0. 
This  Improvement  Is  close  to  fifty  percent  of  the  total 
correlation  energy  of  0,042  a.u.  [5]-   If  we  apply  a 
unitary  transformation  to  get  us  to  the  best  Hartree 
state  then  our  uncorrelated  starting  base  will  be 
-  2.86167  a.u.  [7]  so  our  energy  value  will  be  around 
-2.88167  with  the  use  of  only  two  variable  parameters. 
This  is  probably  the  best  value  that  can  be  obtained  on 
the  basis  of  a  unitary  transformation  approach,  and  only 
two  parameters.   It  can  be  seen  immediately  by  using  a 
natural  spin  orbital  expansion  for  the  first  order  density 
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matrix  I.e. 

that  the  t]    ,    which  are  the  eigenvalues  of  the  \~    matrix 
used  as  an  operator  with  Integration,  cannot  be  altered 
by  a  one  body  unitary  transformation.   This  In  part 
explains  our  Inability  to  obtain  the  total  correlation 
energy.   As  noted  previously  a  comparison  of  the  magnitudes 
of  the  harmonic  oscillator  density  matrix  eigenvalues  with 
those  of  helium,  as  obtained  by  Sobol  [8],  showed  the 
magnitudes  of  the  eigenvalues  In  both  cases  to  be 
fairly  close. 

The  unitary  transformation  used  in  this  calculation 
was  constructed  to  transform  the  nonlnteracting  harmonic 
oscillator  ground  state  into  the  nonlnteracting  helium 
ground  state.   However  it  did  not  carry  out  the  same 
correspondence  to  higher  states.   That  this  fact  did  not 
affect  the  results  much  will  be  seen  by  a  comparison  with 
the  results  of  our  completely  Independent  calculation  on 
helium  using  the  oscillator  based  density  matrix,  with 
only  two  parameters,  which  yields  the  same  minimum  value 
of  E  =  -  2.868  a.u. 

Given  a  good  closed  form  of  U,  calculations  for  large 
atoms  can  be  carried  out  in  the  same  manner.   The  U  used 
in  the  present  helium  calculation  is  not  suitable  for 
larger  systems.   However,  the  feasablllty  and  effective- 
ness of  this  approach  have  been  demonstrated. 


55 


Summary 

We  have  presented  approaches  for  Improving  the 
calculation  of  energy  values  for  atoms  by  accounting 
for  the  inter-particle  correlation.   The  projection  of 
symmetry  breaking  intrinsic  states  was  studied  in  this 
connection.   The  more  powerful  approach  of  intrinsic 
state  density  matrices  was  then  developed.   Oscillator 
based  density  matrices  were  then  studied  and  a  useful 
approximation  technique  was  developed  for  their  applica- 
tion.  The  formalism  was  applied  with  the  use  of  an 
additional  transformation  to  the  case  of  helium  to  show 
the  Improvement  resulting  from  it. 

A  direct  harmonic  atom  approach  was  then  studied. 
Unitary  transformations  allowed  one  to  obtain  a  large 
class  of  one-body  density  matrices,  and  guaranteed 
realizable  two-body  density  matrices.   This  formalism 
was  also  tested  on  helium  and  shown  to  give  results 
similar  to  those  obtained  by  our  previous  technique. 
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Appendix  I.   Details  of  the  Application  to  Helium  of  the 
Ansatz  Wave  Function  Arising  from  Our  Transformation  of 
the  Oscillator-Based  Density  Matrix. 

We  here  present  details  of  the  application  to  helium 
of  the  ansatz  wave  function  arising  from  our  transformation 
of  the  oscillator-based  density  matrix.   Our  trial  wave 
function  Is  given  by 


fl.l)  V/  =  Ne 


with  the  Hamlltonlan 


» 


(1.2)     H  =  -  |^(V>2)  _  |ef  _  |ef  ^  e^ 


2m  12'    r^     r^    ^i  2 


In  the  p  =  0  case,  we  have  the  simple 
(1.5)     ^o=i(^)^-  2   12^ 


which  yields  a  minimum  energy  of  E  =  -2.8476565  a.u.  for 
ctQ  =  3.37500.   We  note  at  this  point  that  the  best  value 
obtained  for  the  solution  of  the  non-relativlstlc 
Schrodlnger  equation  is  E  =  -2.9037  a.u. 

In  carrying  out  our  calculations  the  transformations 
that  It  is  convenient  for  us  to  use  are  as  follows: 
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(1.4)  r  =  p  p    ,     x^  =  ppj_ 

-/r  =   p    ,     y?  r  =   p    ; 


"    y? 


1  2 

—    ,      r  =   p    , 


(1.5) 


^5  2p^    '  ^ 


^X.  pj_pj 


+   5.  .    p 
p  ij    ^ 


These  relations   can  be   checked  using 


(1.6) 


dp  J  "SSTT 
^1       .1 


=  1 


The  Jacoblan  for  transforming  from  the  x  ->  p  system. 


(1.7) 


=  ^(x  y  z)   ^  2  D^ 


^(P^PyPz^ 


(1.8) 

So 

(1.9) 


S  ^  ^-  ^Pl   a 
^     1   ^  ^ 


V^  /(r-Lr2).V^^(r^r2)  d^r-^  d^r2 


:—  .--  ^Pl  ^  S   ,*^  2  2^   a   ,/  2  2, 


•  Ja^b  -^^Pa  -S^Pb  ' 

dp  •  dp  .    5  .  . 

(i.io)  zz^^.  =  -^- 


Pl  =  Pal  • 


3PiP_ 


^k  "^k 


2    hT^ 
P     ^P 
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a/  2^  2,   6-   - 

,^,,,      S     -  2(Pa+Pb)-  2Pa'Pb 

(1. 11)  -T- —  e 

^Pal 

for  our  original  ^(r,rp  )  =  N  e 
for  helium. 

Our  total  expression  for  Equation  (1.9)  now  is 

2   2    — 

2  rr    "^^^Pa'^b^'^Pa'Pb 

(1.12)  ^       J   ^  a   D  . 


[aw4  -  i  ^^4^'  -  -P  ^  1  pIpI  ^\^\ 


a       "^a  "^a 

1 


(The  parts  in  brackets  will  be  considered  separately  and 
designated  as  la,  Id,  lb,  Ic  respectively.) 

Now  consider  the  integrals  needed  to  evaluate  the 
energy .   Let 

^  "  -^^Pa+Pb^'PPa'Pb 

(1.13)  V  ^  ®  ' 

and  2  p 

rr-a(po+Ph)  4  4 

(1.14)  l^ia, a, ^)=JJe  Pa  Pb  ^^"^  ^^  PaPb^^Pa^Pp  ' 

so  our  normalization  condition  is 

(1.15)  1  =  n2  Jk^p  4  p5  p5  dp^  dp^  dO^  dO^ 

p 
=  n2.4.  i^  Ijj(a,a,p)   . 
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In  atomic  units  (e=h=m=l,  aQ=l)  the  energy 
(1.16)     E  =  f  V ^11/*^ -^f  61   ds    -    2  Z    I    f*  -^  f  61   6s    + 


+  I   f*   TT^  ^  dl  ds  , 
^1  2 


the  parts  of  which  we  designate  as  1,  2,  3. 

Part  3  we  consider  separately;  parts  1  and  2  can  be 
written  out  as  follows,  with  dV^^  =  dp^  dp^  dO^^  dCi^, 

la   =  a^/k    ; 

lb   =  -  5  pV  j5rK^p(p-^.p-^)2  p^p5  dV^^  , 

ic   =apN2  ffK^^(p^'P^)    pIp^  ^\^    ; 

Id   =pV  jJyl^^  pIpI^\^  ; 

2    =2Z.8n2JJk^P  P^P^d\b  • 


The  angular  integrations  are  carried  out  using  the 
following  equations: 

(1.17)   //  e"""i''2'  ao^^  ao^^  =  (M^  ^    -. 

1  Sp 

r>  p  -B(e^  'e^)  _  _  h-TT  9 

(1. 18)     Me    -^  "^     ^l^2^^e     ^^e     "  (^)  sinh  B  - 

'   '  -  IM!  co.h  B  : 

a 
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(1.19) 


^1  ^2 


-B(e  -e  )  _   _   p 

Id  ^1   ^2 


-i^^  '^  (-%  +  1)  slnh  B-  I  cosh  B, 


The  resulting  integrals  can  then  be  completed  by  using 
the  results  of  Table  I. 

We  now  consider  the  integration  integral  (3): 

,2 


N 


-aCr^+r^)-^  V^^^   ^l'^2 


^12  d^r-^  d^r^  . 


The  angular  Integrals  will  depend  on  the  sign  of  p 
We  have 


(1.20) 


a 


Dl 


iH. 


dO^  d02 


Sir 


(r^r2) 


'12 

2  ^(r^  +  r2)/2  VrTr^    /-g"    r    +x 

57^^  VlFT  J    "        '^' 


where   c 


Let 


^/ 


^ 


2  /r^r2 


•  |r-,  -rpl    and   d 


^ 


2  /r^r2 


^1+^2 


(1.21)      D(a)    =    /    e 
0 


^2  z?'        2 

^     dx   ;      Ma)    =   /   e  ^     dx   ; 

0 


^(a) 


D(a)   1 
Ma) 
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Use  N^  =  p/4(47r)^Ij,   (for  ^   ^  0)    with  the  scalar  transforma- 

2 
tlon  r  =  p  ,  dr  =  2p  dp  . 


3a) 

b) 


21. 


N 


-a(p^+P^)  3  3  


+lpl(Pa+Pb)/2p3Pt 


^Pa   ^Pb 


BP 


BN 


t£(vK    IPa^bl)    -f  (/S^    "Pa-PbD    1 


2paPb 


a^^b 


For    (3a),    use   the   definitions 


-a^     ^    .2 


^(a)  =  e     Te^   dt  -  DAWF(a);   ERP(a)  =   ^'^     J 
0  0 


2  P   ,2 

e  ^   dt 


2 


.  ^    e^   fl    rr    ^   3,/ r  '^^^Pa+Pb^-^PaPb/^.^p, 

5^)  =T^j2jJ    PaPb>'^IPbf^  '^^^^^ 

-a(p^+p^)-pp^p^ 
-  e  c//(BN)]  dp^  dp^   , 


and 


e^  fK 
^b)  =  —  /^ 
''"  ^    I   V  2 


//Pa\^>^Ie-'^a%^''-^^^^^.H.,BP, 


-  e 


-a(p^+p2)^p^p^ 


ERF(BN)]  dp^  dp^  . 
These  expressions  were  evaluated  numerically, 
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-  i(a^x^+b^y^-cabxy) 

Let   K^bc  =   e      ^ 

with  a>0,   b>0,    |c|    <1. 
oo    oo 

0   0 


Table  I 


cos        c 


0      -        1  r        1  C    COS         C    -1 

(ab)  1-c  (1-c    ]    ' 


-1  ,2      „_-l 


/9      _          1          r                 3C  .     cos          C             3c       COS          c    1 

/o             1        r          4          ,  15c              9c   cos        c  _  3-5c^cos     c-, 

/9             1        r              55c  7-15  c^       9  cos"-*-   c 

P        -1  4-1 

.    2 •  5 •  9   c   cos        c  ,    3-5-7   c      cos        c    -, 
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Appendix  II.  Details  of  the  Application  to  Helium  of 
the  Harmonic  Atom  Approach. 

We  here  present  details  of  the  application  to  helium 
of  the  harmonic  atom  approach.   The  expressions  In  Eq.  (6.22) 
are  expanded  by  the  use  of  (6.26).   Using  considerations 
of  symmetry  In  the  hamlltonlan  the  following  breakup  Is 
obtained: 

r^*(l"2")  [  Tu  U  H  U  U  dl  d2]  ^(1'2')  dl'd2'dl"d2" 

r/(l"2")  H(l"2")  ^(1"2")  dl"  d2"  I 

-  2  r^*(l  2')  H(l)  A(l  1')  ^(1'2')  dl  dl'  d2  •  Ila 

-  2  r^*(l"2')  A(1"1)H(1)^(1  2')  dl  dl"  d2'  lib 

+  2  r^*(l"2')  A(l"l')  H(l)  ^(1'2')  dl"  dl'  d2'  lie 

_  4e^  f  f*{l   2)  A(2  2')  -^  ^(1  2')  dl  d2  d2 '  III 

J  ^12 

+  2e^  fi^*{l   2")  A(2"2)  A(2  2')  -^  ^(1  2')  dl  d2  d2 '  d2"   IV 
J  ^12 

+  2e^  r^*(l  2)  A(l  1')  A(2  2')  -^  ^(1'2')  dl  dl '  d2  d2 '   V  LI 
J  ^12 

+  2e^  f  i/*{l   2")  A(l  1')  A(2"2)  -^  ^(1'2)  dl  dl '  d2  d2"   V  L2 
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ke'^    [i/*{l   2")  A(l'l)  A(2"2)  A(2  2')  -^   ^(1'2') 


12 
dl  dl'  d2  d2'  d2"        V  L5 


+  e^  r/(l"2")  A(l"l)  A(l  1')  -^  A(2"2)  A(2  2 


12 
V'(1'2')  dl  dl'  dl"  d2  d2'  d2"    V  L4 


where 


TT 

2 

(II. 2)     H(l  2)  =  H(l)  +  H(2)  +  ^ 


^1  2 


2 

(11.3)  H(l)  =   -^  civl)    -   z  e^  S(^)  ; 

(11. 4)  ,5)  ^  j^^4_(l)0/L(l)  dl  ; 

(11. 5)  A(l  2)  =  A(l)  A(2)   . 

Some  parts  of  these  Integrals  are  not  integrable 
analytically.   So  the  integrals  will  be  reduced  to  double 
Integrals  for  which  I  have  developed  Integration  routines 
for  the  computer  on  the  basis  of  a  double  Simpsons  rule. 
A  part  of  this  reduction  Is  most  easily  evident  by  examining 
(6.23). 
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Let 

(11. 6)  B(l)  =/A(2)  ^(1  2)  d2 

(11. 7)  H^    2)  =  Tud  1')  U(2  2')  ^(1'2')  dl'  d2' 

=  ^(1  2)  -  A  (2)3(1)  -A(i)B(2) 

+  A(l  2)  /  /V(2'  )  3(2'  )  d2'  . 

^^ . ' 

TAB 
Now  take  ^{1   2)  •  (^(1  2).   When  this  expression  is  Integrated 
or  Is  multiplied  by  a  symmetric  term  like  l/r.,  „  and 
then  Integrated  we  can  by  symmetry  reduce  the  Integrals  to 
sevel  terms.   These  seven  terms  can  be  Identified  with  all 
the  l/r-,  p  terms  of  our  previous  breakup  and  could  be 
obtained  by  further  integration. 

(11. 8)  ^(1  2)  ij{l    2)  -  4^(1  2)/4(l)B(2)  +  2-^(1  2)A(1  2)TAB 
[Part  of  I]        [III]  [V  LI] 

+  2[A(2)&(1)]^  +  2A(2)B(1)A(1)B(2)-^A(2)B(1)A(12)TAB 
[IV]  [V  L2]  [V  L3] 

+  A(l  2)A(1  2)(TAB)^    ^[V  L4]  . 

We  will  consider  the  reduction  of  the  angular  integrals. 

For  dealing  with  most  of  the  kinetic  energy  (II)  terms, 

it  is  more  convenient  to  work  with  the  breakup  of  Eq.  (6.22). 
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Note  in  the  following  review  of  II  terms,  we 
Interchange  orders  of  integration.   It  has  been  verified 
that  this  is  proper. 

(11. 9)  r^*{l    2)  H(l  2)  (t)(l  2)  dl  d2 

=  /  ^*(1"2")  U*(l"l)  U*(2"2)H(1  2)U(1  l')U(2  2')^(l'2') 

dl  d2  dl'  d2'  dl"  d2"  . 
For  the  single  particle  part,  i.e.  H(l),  we  can  use 

(11.10)  j\J*{2"2)    U(2  2')  d2  =  5(2"-2')  , 
to  simplify  integration.   Now  in 

(11.11)  r^*(l"2")  U*(l"l)  H(l)  U(l  1')  i//(l'2")  dl  dl'  dl"d2' 
with  ^(1  2)  =  N  exp  [-  ^(r^+r^)  +  Br-j_.r2^'  '^ake  rd2", 

(11.12)  r/(l"2")  ^(1'2")  d2"  = 

J  2  2 

(11.13)  UHU   =   5H5-5HA-AH54-AHA. 
Consider  the  last  three  terms  a',  b',  c'. 
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Take  dl" 

2  2 

2  ^  3/2p  -  -|(A-  ^)(r  +r  )+  1^  r  .? 
(a')  =  -N  (|)    /  e       2A    i   i     2A   i   1 


H(l)  A(l  1 ' )  dl  dl 


(Note  H  act   on  quantities  to  the  right  only.) 


(b-)  = 


2  TT  5/2  r 
N^(J)   J   A(l"l)  H(l) 

-  |(A-  ||)(r^„+r2)+  g^  V 


(C) 


2  TT  5/2  r 

+  n2(|)    /  A(l"l')  e 


2A  '1"  "1 

1,„   b2\  ,  2  ^  2  .^  B^  - 

2^    2A^^^1"^^1'^   2A  ^1"'^1' 


dl  dl"  ; 


dl"  dl'  yy^(l)H(l)^(l)  dl 


Use 


(11.14)   V^  e    ^    ^   ^ 


-6a+4a  r-,  +  p.  rp  -  ^apr.,  -rp, 


L 


•  e 


-ar-j^+  pr-j^  -r^ 


a 


5/^   2__2  ^  ^  _   2  ^1 


(11.15)  yi-(uv)  V^^(l)  =  (^)    (a2r^--a)e 


a-  r^^^^^^r2b   ,  2^ 
+  ( — )    ( —  -be 
1^  r-. 


-br 


Note  the  6H5  contribution  coincides  with  Eq .  (6.l6),  which 
has  beeri  calculated  previously  and  is  our  entire  I  term. 
We  now  define  a  number  of  quantities  and  present  some 
integrals  to  be  used  in  reducing  the  order  of  the  II  integral 
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Let 

2 

(11.16)  EXB  =  e 

1  , ,  2,  2>   B^,      .2 
-   Ts  A(r  +r  )+  -[TT-Cr  -r  ) 

(11.17)  EXM  =  e   ^     ^   ^    4A   1   2 


(11.18)       DLA  =  r^V^  4  (1) 


(11.19)  -  TCB  =  H  =  /  A  (1)  H(l)  A  (1)  dl 

=  ^/[-2   A(l)  r^  DLA  -  8r^  A^CD]- 
The  resulting  sum  of  the  II  terms  a,b  c,  is 

(11.20)   2|_^_  f   /^   (2)  r J(EXB-EXM)(0.5  *  DLA  +  A  (D 

2  4 

B  A,  r,  r-//\2  T^2  ,   B  >,_;^ 


•  [4.-r^(TCB-0.5  ^  A  +  ^)+  0.5(A^-B^+  -^)r^ 

+  ^r^rj)  -  (0.5  B^-  ^)r^r„(EXB +EXM)  >dr,  dr„ 
8A    "^  ^  ^A"^   ^  ^  r   1    ^ 

For  B  =  0  limit 

1     2    2 
(11.21)       M-EXM^   ^-2^(^1-^-2^,1^, 

B^  ^ 

(Note  the  above  results  are  in  my  "bar"  notation  where 
N^  =  ir^N^.  A  =  TT  \  >    DLA  =  ir   DLA.   These  definitions  can 


used  to  convert  to  ordinary  notation.) 

The  following  equations  have  been  used: 

—   —  P 

P    pr  -r  ihir)      slnh  p  r-,r2 

(11.22)  /  e   ^   ^  dO.  dOp  =  

(11.23)  /?l-^2^   ^   ^  dO^  dOg   =^-(11.22) 

.   2   slnh  p  r  r     ,i^    .2 
==  -  (-11)"^  ^^  +  ^^J^    cosh  p  r,r„ 

00       p  2 

(11.24)  /  X  e-^^  slnh  x  R  dx  =  ^^^  /|  R  e^  /^^ 

0 

f   -Ar^  ^_   _   B?  •?  ,+  C?2-^l 

(11.25)  j  e   ^  r^  6^r^J   v^-r^   e   ^   ^  dO. 

0  r 


j  e     r2  dr2  ^  j  e         d02 
0  r 


with  R  =  |Br-^,+  Cr-^l  =  ^(Br^  +  Cr^)-(Br^,+  Cr.^)  ;  we 
apply  (11.22)  and  assume  dr2  Integration  can  be  carried 
out  before  taking 


can 


R^  ./o        o.n        2 


;,            3/2      -  ^          ^3/2      2(C  rf'+B  r    -r    ,) 
(11.26)      ^    (a)  e  =    (a^  JfA 

(Br^,4€?^) •(Br^,+Cr^)/4A 
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^  -Ar   1,      n    Br  e  -e 
(11.27)   /  e   ""  r^  dr2  /  e   "^  "^   60^  = 


S    ^  3/2   B^AA 


Note  for  B  =  C  in  (11.26),  let  B  =  C  in  result  and  get 
(P     2^  (^1  -^^1-^1')  ^ 


Our  II  terms  now  give 
Ila  =  -  2  r/(l  2')A(1')  H(1)A(1)  ^(1'2')  dl  dl  •  d2 '  . 

^  _  _|2    N     r  (gxB  -  EXlVDr^r,  ,  A  (!')  H(l)  Ad)  dr..dr 
TT^/^  -/A  B^  ^  ^  ^ 

=  +-22  _^T     r(EXB  -  EXM)r   AH^)- 
TT^/^  Va  B^  J  ^ 

.[|  r^  V^  +  2]   A  (1)  dr-^  dr^ .  • 
lib  =  -  2  r/,(l'2')  A(l'l)  H(1)  ^(1  2')  dl  dl'  d2' 
-^;;372  ^32   2j  A(l  Dr^r^e 
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•  I=i"h  Is  r,r^,^^  -  I  +  ^  +  |(A^-B^  ^)r^+  -S!  r^ , 


4A 


4         2 

g)  cosh  2^  r^r-^  ,1  dr^  dr^  , 


lie  =   +  2   E        Tp    (1'2')  A(l  1')  ^(1  2')  dl  dl'  d2 ' 


-  """372 


N 


TT' 


Va  b' 


(EXB-EXM)  r^r-^  A   (1  1')  dr^  dr^  , 


H. 


We  now  must  evaluate  ]3(1)  =  /A(2)^(1  2)  d2,  which  is 
used  in  the  reduction  of  the  l/r-,  „  terms.   Using  our  table 
of  Integrals,  we  get 


A   2 
2  ^2 


(11.28)   B(2)  =  N 


A   2^  ^-   - 
-2  ^i+Br^-r^ 


TT  Y  TT 


-br-, 
e    -^>  dl  ; 

J 


P    Br-,  -rp        J. 
(11.29)   /  e     ^  dl  =  ^^   slnh  (Br^r^)  ; 


Br^rg 


(11.50)   B(2) 


N  e 


-A  ^2 
2   2 


1-    VqUq  d2' 


„  3A  o^  3/2 


IT 


'A4a 


"  2  ^A  ^    VV   rF 


^    [(b+Br^' 
TT   Brp       2' 
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(b+Br  )  /2A       b+Br   -        (b-Br^)V2A 

•  e  erfc  ( ^)-(b-Br^)e 

/2A        ^ 

b-Br^   1 

•  erfc  ( -)]  y   . 

/2A   J 

The  complementary  error  function  is  denoted  by  erfc(x). 
Let-J?(BR)  =—[...].   At  Br  =  0,  we  must  take  a  Taylor 
expansion  or  do  the  integration  for  this  point  specifically 
The  result  obtained  either  way  is 

2   b  /2A  

^(0)  =  2  [(1  +^)  e      erfc(^^)  -  /^  b]  . 

Note  in  the  B  =  0  limit  when  i/{l   2)  =  Uq(1)uq(2),  using 
^{0) ,    and  rearranging  constants,  this  does  reduce  to 
the  proper  limit  of   13(2)  =  Uq(2)  ^1-  /'  v^u^  d2'  . 

We  now  can  give  a  review  of  the  reduced  form  of 
the  rest  of  the  integrals  that  were  used. 

f/il   2)   (2)  -^—  B(l)  dl  d2  111-4 

J  ^12 

2  [     (1)  A(2  2)/r^  2  ^d)  dl  d2  IV 


2        ip    (1   2)    A(l  2)/r^2'^^B  ^^  ^2  , 

TAB  =  rA(l'2')^(l'2)  dl'd2  V  LI 

2  r  5  (1)  A(l  2)/r^  2  '-^(2)  dl  d2  V  L2 
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2     r     A(ll)  A(22)     ^^    ^2  V  Lif 


(TAB) 


For   further  reduction,    use 


P    do,    dO„  o^2 


A,    2,    2^ 
(II. 32)       /   ^    (1   2)/r-|^2     dO^   dOg  =   ^TT  N   e      ^      ^      ""    • 

•f^  (|)'^'t  1(71  [-1-2^ )-i(/l  1-1-2')^ 

|(r^+r2) 
•    e 

For  B  small  or  for  negative  B,  we  replace 

Jim  A  i''i^'^2i'  -i</l  1^^1-21']  ■ 

where  J  is  the  error  function  by 

2   CD   (-B)^[(r^+r2)2"+l-lr^-r2|2^+^] 


yi   0  2^  n:(2n+l) 
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